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Abstract Recent experimental results on the static or 
quasistatic response of granular materials have been in- 
terpreted to suggest the inapplicability of the traditional 
engineering approaches, which are based on elasto-plastic 
models (which are elliptic in nature). Propagating (hyper- 
bolic) or diffusive (parabolic) models have been proposed 
to replace the 'old' models. Since several recent experi- 
ments were performed on small systems, one should not 
really be surprised that (continuum) elasticity, a macro- 
scopic theory, is not directly applicable, and should be re- 
placed by a grain-scale ("microscopic") description. Such 
a description concerns the interparticle forces, while a 
macroscopic description is given in terms of the stress field. 
These descriptions are related, but not equivalent, and 
the distinction is important in interpreting the experimen- 
tal results. There are indications that at least some large 
scale properties of granular assemblies can be described 
by elasticity, although not necessarily its isotropic version. 
The purely repulsive interparticle forces (in non-cohesive 
materials) may lead to modifications of the contact net- 
work upon the application of external forces, which may 
strongly affect the anisotropy of the system. This effect is 
expected to be small (in non-isostatic systems) for small 
applied forces and for pre-stressed systems (in particular 
for disordered systems). Otherwise, it may be accounted 
for using a nonlinear, incrementally elastic model, with 
stress-history dependent elastic moduli. Although many 
features of the experiments may be reproduced using mod- 
els of frictionless particles, results demonstrating the im- 
portance of accounting for friction are presented. 

Received: February 2, 2008 

Chay Goldenberg^ and Isaac Goldhirsch'^ 

^School of Physics and Astronomy 
Tel- Aviv University 
Ramat-Aviv, Tel-Aviv 69978, Israel 
e-mail: chayg@post.tau.ac.il 

■^Department of Fluid Mechanics and Heat Transfer 

Faculty of Engineering 

Tel- Aviv University 

Ramat-Aviv, Tel-Aviv 69978, Israel 

e-mail: isaacOeng. tau. ac . il 



Key words. Granular response - elasticity - coarse-graining 



We appreciate helpful interactions with R. P. Behringer, 
J. Geng, E. Clement, D. Serero, H. Jaeger, T. A. Witten, N. 
Mueggenburg and J.-N. Roux. Support from the U.S.-Israel 
Binational Science Foundation (BSF), INTAS and the Israel 
Science Foundation (ISF) is gratefully acknowledged. 



Introduction 

The modeling of granular materials has been a subject 
of ongoing research in the engineering community (see 
e.g., [1]). In recent years, this subject has found renewed 
interest among physicists [2-5] (having been studied in the 
distant past by great physicists such as Coulomb, Faraday, 
Reynolds and others). 

The behavior of "granular gases" , which are obtained 
by e.g., sufficiently strong shaking or shearing (so that 
the material behavior is dominated by interparticle col- 
lisions), has been quite successfully modeled using ap- 
proaches based on extensions of the kinetic theory of gases 
[6]. However, the behavior of dense granular matter, which 
is dominated by prolonged interparticle contact, has proven 
more difficult for modeling. For the description of the 
quasi-static behavior, elasto-plastic models are commonly 
used by engineers [7, 8]. 

This paper is concerned with the static behavior of 
granular systems. In elasto-plastic models, one often uses 
(linear) elasticity below yield (although parts of a static 
system are sometime assumed to be at incipient yield [8] ) . 
However, in recent years a very different class of models 
has been proposed for describing the statics of granular 
materials, based on the notion of "force propagation" , sug- 
gested by the observation of force chains in experiments 
on granular materials [9], as well as simulations [10,11]. 
These models (see e.g., [12-15]) typically yield hyperbohc 
partial differential equations for the stress field, in contrast 
with the elliptic, non-propagating nature of the classical 
equations of static elasticity. It has been claimed that the 
hyperbolic description tends to an elastic-like one at large 
scales [16-18] (however, the physical interpretation of the 
macroscopic fields in this case is not clear). 

Recently, the response of granular slabs resting on a 
horizontal floor to a 'point force' applied at the center of 
the top of the system has been studied experimentally [19- 
25]. In [19,22,23], the intergrain force distribution has 
been measured in two-dimensional (2D) systems as a func- 
tion of vertical and horizontal distance from the point of 
application of the force. In [25], the particle displacements 
for similar 2D systems have been measured. In [20, 21, 24], 
the vertical force acting on the floor has been measured in 
three-dimensional (3D) systems. Prominent force chains 
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have been observed in ordered 2D systems; these force 
chains fade out with increasing disorder. For pentagonal 
particles in 2D arrangements the measured force distribu- 
tion is single peaked and the width of the peak is linearly 
related to the vertical distance, in conformity with elas- 
ticity. The results for cuboidal particles obtained in [19] 
appear to suggest a parabolic behavior, consistent with a 
diffusive model, although the systems studied were quite 
small. In [25], the width of the measured distribution of 
displacements, as a function of the vertical distance from 
the particle which is directly displaced, follows a square 
root dependence (as expected from a diffusive model) for 
small distances of a few particle diameters, crossing over 
to a linear dependence at larger distances (consistent with 
an elastic description). Ordered 3D packings exhibit mul- 
tiple force peaks for shallow systems [24] and less structure 
for deeper ones. Somewhat larger (in terms of number of 
particles), disordered 3D systems [20,21] exhibit a single 
peak in the force distribution measured at the floor, whose 
width is proportional to the depth of the system. 

The experimental evidence appears to be contradic- 
tory: different experiments seem to support fundamentally 
different descriptions of the response of granular mate- 
rials (in the case of 2D systems, it has been suggested 
that there may be a crossover from a hyperbolic to an el- 
liptic behavior with increasing disorder [22]). The thesis 
presented in this paper is that these seemingly contradic- 
tory experimental results (and theoretical explanations) 
are not necessarily at odds with each other. This the- 
sis is based on the observation that most of the studies 
(perhaps all) rejecting the elliptic description have been 
devoted to small systems, of the size of a few dozen of par- 
ticle diameters at most, whereas many engineering stud- 
ies consider rather large granular systems. Since elasticity 
and other macroscopic descriptions are not valid on small 
scales, at which local anisotropics and randomness play 
a major role, one should not be surprised that such de- 
scriptions fail on small scales. Indeed, simulations [26, 27] 
reveal the existence of a crossover from microscopic to 
macroscopic behavior of granular assemblies (as well as 
other systems [28]) as a function of system size or resolu- 
tion. We argue that such a crossover is observed in some 
of the experiments mentioned above. Strictly isostatic sys- 
tems [29] have been shown to be described by hyperbolic 
stress equations [15,30], and numerical simulations sug- 
gest that systems of frictionless spherical particles ap- 
proach isostaticity in the limit of infinite rigidity [31]). 
However, we argue that since real granular systems have fi- 
nite rigidity and usually experience frictional interactions, 
they cannot be generically isostatic (the same presiimably 
holds even for frictionless non-spherical grains). The iso- 
static limit is a singular case, whose physical consequences 
for real systems are at best unclear. Therefore the contro- 
versy surrounding the correct description of granular stat- 
ics is mostly a question concerning the behavior of small 
granular systems. The latter require a grain-scale ("mi- 
croscopic" ) description, rather than a macroscopic one. 

A second point stressed below is the distinction be- 
tween force and stress. Whereas interparticle forces can 
exhibit force chains which look like they contradict elas- 
ticity, the latter does not describe the nature of the forces 



but rather that of the stress field. The stress field involves 
an averaging over the forces (whose result is resolution 
dependent) and leads to less pronounced structure than 
the underlying force field. The small scale structure of the 
interparticle forces cannot be taken to consist an argu- 
ment against an elliptic description or in favor of it, since 
it relates to small scales and it does not deal with the 
objects with which elasticity or plasticity are concerned. 
The large scale response of granular packing is shown to be 
consistent with a (possibly anisotropic) elastic description. 
The fact that in non-cohesive granular materials there is 
no significant attraction among the particles may lead to 
modifications of the contact network, which may strongly 
affect the anisotropy of the system. This effect is expected 
to be small for small applied forces (for non-isostatic sys- 
tems) and for pre-stressed systems, in particular for disor- 
dered systems. Otherwise, it may be accounted for using a 
nonlinear, incrementally elastic model, with stress-history 
dependent elastic moduli. 

The third point made in this paper is that while models 
employing frictionless particles can reproduce some prop- 
erties of granular packings, friction can be of utmost im- 
portance for the description of granular matter (a rather 
intuitive fact). Results demonstrating the importance of 
accounting for frictional interactions are presented in Sec. 5 
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The microscopic picture: forces 

In attempting to describe granular materials in terms of 
continuum mechanics, by analogy to "regular", atomic 
materials, one usually considers the "microscopic" scale 
to be that of the individual particles (whose internal dy- 
namics should be well des(;ribcd by continuum mechanics) . 

One of the simplest granular systems is a collection of 
frictionless spherical particles. A typical microscopic (par- 
ticle scale) description of such a system is given by the 
particle's radii, {Ri}, their masses, {nii}, center of mass 
positions, {ri{t)}, and velocities, {vi{t)}, at time t. It is 
typically assumed (e.g., in the context of simulations of 
granular materials [10, 32, 33]) that the particles are quite 
rigid, so that the interaction between two particles (in 
the frictionless case) depends only on their respective dis- 
tance, or, more conveniently, on their imaginary overlap 
C„(t) = Ri+Rj-\r^j{t)\, where (t) ee nit) - rjit). The 
contact interactions are usually modeled by treating the 
particles as macroscopic objcicts, dc!S(;ribed by the equa- 
tions of continuum mechanics (see e.g., [34,35]). For two 
frictionless elastic spheres, a classical result by Hertz (see 
e.g., [36]) is that the force is proportional to S,"^^^, while 
for cylinders, it is linear in the overlap. For noncohesive 
particles, only repulsive forces are possible. Even for fric- 
tionless particles, internal dissipation as in e.g., viscoelas- 
tic particles, gives rise to a dependence of the force on 
the relative velocity ^ as well (for some examples of force 
schemes commonly used in simulations, see e.g., [37-39]). 

The interparticle forces for a given configuration of 
such particles subject to given boundary conditions (e.g., 
specified displacements of the particles on the boundary, 
or forces applied to them) and body forces such as gravity 
can be determined, for a static system, using the equations 



3 



of equilibrium (Newton's laws) and the force-displacement 
relation. We reiterate that when full force laws for par- 
ticle interactions are known or modeled, the statics and 
dynamics of the system arc fully determined (they may 
be history-dependent for history-dependent force laws, as 
commonly used for frictional interactions). 

In the case of frictionless isostatic systems (in which 
the mean coordination number is exactly z = 2d, where d 
is the dimension of the system) the forces can be deter- 
mined from the equations of equilibrium alone (and are 
therefore independent of the force-displacement law; how- 
ever, the particle displacements certainly depend on this 
law). It has been suggested [40,41] that frictionless granu- 
lar systems become isostatic in the limit of infinite rigidity 
(giving rise to a macroscopic behavior which is very differ- 
ent from elasticity), and this appears to be borne out by 
numerical simulations [31,42]. However, the relevance of 
this limit to real materials is questionable, since real ma- 
terials cannot be infinitely rigid. Any additional contacts 
created if the rigidity is allowed to be finite will render 
the system hyperstatic (so that there is a "phase transi- 
tion" to an isostatic behavior only at infinite rigidity [40]). 
The rigidity should of course be compared to the confining 
forces or body forces (in a system under gravity and con- 
fined by walls, the confining force is related to gravity). 
If the confining forces are very small, the system would 
indeed be expected to be close to marginal stability. As 
mentioned above, the static indeterminacy associated with 
hyperstatic systems simply means that the equations of 
equilibrium are insufficient for determining the forces, so 
that additional equations (e.g., force- displacement laws) 
are required. Static indeterminacy does not mean that 
there's no unique solution for the forces in a real system. 
A similar situation occurs on the macroscopic, continuum 
level (see e.g., [8]). The rigid limit can be approached in 
many different ways (e.g., the stiffness of each interparti- 
cle contact may be different) , and, even if assuming that 
the same (isostatic) contact network is obtained for differ- 
ent distribution of the interparticle stiffness, yielding the 
same interparticle forces, the particle displacements will 
certainly be different, hence the rigid limit in not unique, 
at least in this sense. 

In several experiments, photoelastic particles were used 
in order to measure the stress in granular systems [9, 
19,22,23,43]. These measurements probe the intraparti- 
cle stress, i.e., the stress inside each particle. Following 
the above, these should be interpreted as measurements 
of microscopic fields (the macroscopic description of gran- 
ular systems regards the particles as microscopic, and does 
not resolve any details below the particle scale). The mi- 
croscopic fields corresponding to these measurements are 
the interparticle forces, which can be deduced from these 
internal stress measurements (as described in [23]). As 
mentioned, these forces should be distinguished from the 
"macroscopic" stress field in the system. 

The distribution of force magnitudes in a static gran- 
ular packings is a microscopic quantity which has been 
extensively studied in experiments [44, 45] and simula- 
tions [46]. An exponential behavior of the distribution 
at large forces appears to be quite universal in exper- 
iments on granular systems, independent of the degree 



of disorder [45], the friction coefficient [45], or the rigid- 
ity of the particles [47], and has also been observed in 
simulations of granular systems with different models for 
the interparticle forces (e.g., [48,49]). The universality of 
the force distribution appears to extend to other systems 
such as foams, glasses, colloids etc. (see [50] and refer- 
ences therein) . The exponential tail of the distribution is 
reproduced in simple models such as the (parabolic) q- 
model [51,52]. The distribution for smaller forces appears 
to be less universal, and it has been suggested that the 
appearance of a peak in the force distribution near the 
mean force may signal the onset of jamming or a glass 
transition [53]. 

Interestingly, a qualitatively similar force distribution 
is obtained in purely harmonic networks: Fig. 1 shows the 
force distribution obtained for an ensemble of random net- 
works of linear springs constructed as follows. Points are 
placed on a 2D triangular lattice with spacing d (with 
square-shaped boundaries), and then their x and y coor- 
dinates are randomly displaced by ±0.04cZ. Points whose 
distance is less than 1.02d are connected by linear springs 
(whose equilibrium length is equal to this distance) with 
equal spring constants (this results in an average dilution 
of about 12% of the springs compared to the perfect lat- 
tice). A uniform isotropic compression of 1% is applied 
to the boimdary particles, and the interparticle forces arc 
calculated. The force distribution presented in Fig. 1 is 
obtained from an average over the force histograms of 100 
systems of 1085 particles. The force was normalized by the 
mean force in the ensemble (a similar distribution is ob- 
tained for a normalization by the mean force for c^acli sys- 
tem; the variation in mean force among different systems is 
relatively small, which may indicate that the system is far 
from 'jamming' [53]). The tail of the logarithm of the dis- 
tribution is fit quite well with a line of slope —3.8, similar 
to the slope obtained in experiments on highly compressed 
disordered packings of soft rubber spheres [47] (similar dis- 
tributions were obtained for a scalar harmonic network of 
unequal springs in [48]). For the case of networks with no 
force dilution (the same connectivity as in the perfect lat- 
tice), the force distribution is Gaussian with a half- width 
of a few percent of the mean, i.e., a much narrower distri- 
bution) . These results indicate that a random connectivity 
should be consequential for the force distribution, which 
may be the reason that even for highly compressed disor- 
dered spheres (whose contact network is still disordered), 
the distribution is qualitatively similar to that observed in 
less compressed systems [47] . A similar effect has been ob- 
served in simulations of granular systems under different 
applied pressures [48]. 

The forces in one of the realizations of the ensemble are 
shown in Fig. 2. Force chains are clearly observed (note 
that there are very few tcnsional forces, so that they do 
not significantly affect the force distribution in this case). 
Similar force chains have been observed in a polydisperse 
Lenard-Jones system [28] (incidentally, the concept of a 
force chain is not well-defined: in the case of a homoge- 
neous strain applied to a uniform lattice, the forces are 
equal, so that it is reasonable to define the force chains to 
contain forces whose magnitude is larger than a uniform 
cutoff, e.g., the mean force, as used in Fig. 2; However, 
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Fig. 1. The distribution of force magnitudes in boud-dilutcd 
distorted triangular networks of linear springs (see text). 



for a non-uniformly strained system, e.g., systems sub- 
ject to gravity, in which the mean force increases with 
depth, such a global cutoff makes little sense). The re- 
sults shown in Fig. 2 indicate that force chains are not 
specific to granular systems. Force chains are microscopic 
features of microscopically disordered systems (or even 
inhomogeneously strained ordered systems, as described 
below), and their presence does not necessarily indicate 
any macroscopic inhomogeneity, or inconsistency with a 
macroscopic elliptic, or elastic, description. It is quite cer- 
tain that if one could observe the individual intcrparticle 
forces in atomic systems (which may not be quite well de- 
fined, since a quantum description is appropriate for such 
systems), one would also observe force chains. 

It is important to note that a significant portion of 
the stress (even in a homogeneously strained system) is 
carried by forces which do not belong to the force chains. 
An example is provided by a system of frictionless polydis- 
pcrsc disks (with radii uniformly distributed within 10% of 
the maximum radius) which is confined by side walls and 
a floor, with a uniform force applied to the particles of 
the "top" layer (without gravity). The intcrparticle forces 
are taken to be linear in the overlaps. Fig. 3 shows the 
forces in the system. Fig. 4 shows the fraction of the ap- 
plied vertical force carried by the forces whose magnitude 
is greater than the mean (i.e., those belonging to force 
chains, using the definition mentioned above), compared 
to that carried by all the forces (which is of course equal 
to 1), for forces in horizontal "slices" of the system, as a 
function of the vertical coordinate, z. As seen in Fig. 4, 
only about 80% of the applied force is carried by the force 
chains. Furthermore, the; force carried by the chains fluc- 
tuates with depth, so that the forces in the chains do not 
obey the conditions of force equilibrium. It is therefore 
questionable whether a model which describes the stress 
exclusively in terms of the force chains is justifiable. 

The (near-) universality of the force distribution, in 
particular the fact that it is observed in simulations of ran- 



Fig. 2. The forces in a bond-diluted distorted triangular net- 
work of linear springs. Forces with magnitude larger than the 
mean are indicated by solid lines whose width is proportional 
to the force magnitude; smaller forces axe indicated by thin 
dotted lines. Compressive forces are indicated by gray lines; 
tensile forces by black lines. 



dom systems with harmonic interactions, does not make 
possible the differentiation between different models on 
the basis of the force distribution (in particular, the ob- 
servation of such a distribution does not preclude an el- 
liptic description). The same statement applies to the ob- 
servation of force chains. A more sensitive and direct test 
should be rendered by the response of a granular system 
to inhomogeneous external forcing, such as that provided 
by localized forces. The latter seem to be consistent with 
elasticity, as described below. 



Macroscopic fields and continuum equations in terms 
of microscopic quantities 

Continuum descriptions of materials are often based on 
phenomenological arguments (usually motivated by exper- 
imental findings), rather than on derivations from the un- 
derlying microscopic dynamics. A unique feature of gran- 
ular materials is that due to the typically large sizes of 
the constituents, it is relatively easy to access the "micro- 
scopic" scales experimentally. On the other hand, in most 
practical applications, the number of particles is such that 
a detailed particle-level description becomes intractable, 
and a continuum description is required. The fact that 
experiments on granular systems can yield both micro- 
scopic information and macroscopic information (possibly 
even in the same experiment) is useful to the elucidation 
of the connection between these two descriptions. 

In order to obtain a macroscopic description of a sys- 
tem in terms of the microscopic fields, we employ a spatial 
coarse-graining approach [27,54]. We stress that the only 



Fig. 3. The forces in a system of polydisperse frictionless dislcs with a uniform force applied to the top layer (no gravity). Line 
widths are proportional to the force magnitudes. Left: all forces, right: only the forces whose magnitude is larger than the mean. 



0.8 -^^aA^a-^aa^aA aaaAA 



f/f 



0.6 

app 
0.4 



0.2 



• All forces 

A Forces in chains 



z/d 



10 



15 



Fig. 4. The fraction of the applied vertical force carried by all 
the forces and by the forces whose magnitude is larger than the 
mean (i.e., those belonging to force chains), as a function of the 
vertical coordinate, 2, scaled by the mean particle diameter, d, 
for the system shown in Fig. 3. 



averaging considered here is spatial (the approach can be 
extended to include temporal coarse-graining as well [54] , 
but here we consider static configurations). Since static 
granular packings are typically found in metastable states, 
far from equilibrium, and thermal energy scales are negli- 
gible, such systems do not explore any phase space so that 
it is hard to justify the kind of ensemble average commonly 
used in statistical mechanics. An average over configura- 
tions (i.e., average over different disordered systems which 
are presumed to be prepared in the same way) is com- 
monly performed when analyzing experimental data, due 
to the large fluctuations obtained in many experiments. 
However it is not clear a-priori if self-averaging occurs, 
i.e., that at least for large enough scales the macroscopic 
behavior of a single "typical" realization is the same as 
that of the average behavior over many realizations. Self- 
averaging may be valid for some quantities and not for 
others. Therefore we choose not to assume a-priory that 
any ensemble averaging is justified; instead we relate the 
macroscopic and microscopic fields in a way that is rele- 
vant for single realizations. 



Following [54], define the coarse-grained (CG) mass 
density p{r,t) and momentum density p{r,t) at position 
r and time t as 



P{r,t) = ^mi(t)[r - ri{t)], 

i 

Pir,t) = ^miV,{t)(t)[r ~ r,(i)]. 



(1) 
(2) 



where <^(i2) is a non-negative coarse-graining function (with 
a single maximum at i? = 0) of width w, the coarse- 
graining scale, and / (j){R)dR — 1. 

Upon taking the time derivative of the macroscopic 
fields p and p, performing straightforward algebraic ma- 
nipulations [54] and using Newton's laws, one obtains the 
equation of continuity and the momentum conservation 
equation, respectively: 



P 



Pa 



-div(py) 
d 

drp 



(3) 



E 



fiVaVf, 



where the velocity field is defined by V = p/ p, Greek in- 
dices denote Cartesian coordinates, and the explicit de- 
pendence of the CG fields on r and t has been omitted 
for compactness. Since this paper focuses on the stress 
field, we have omitted the energy equation, which can be 
derived in a similar way [27]. 

In addition to obtaining the standard equations of con- 
tinuum mechanics from microscopic consideration, this 
coarse graining procedure provides an expression for the 
stress tensor aaf3 in terms of the microscopic entities: 



m, v[^{r, t) v[p{r, t)(j){r - r,{t)) 



(4) 



ds(f)[r - ri{t) + sTijit)], 



where v'^{r, t) = Vi{t) — v{r,t) is the fluctuating velocity, 
fij{t) is the force exerted on particle i by particle j, and 

The flrst term in Eq. (4) is the kinetic stress (which 
vanishes for static configurations), and the second term is 
known as the contact stress. Note that the standard Born- 
Huang expression [55]: CTq/j = -w^ijev-i^j hanjis is 
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equivalent to the expression for the contact stress in Eq. (4) 
if the coarse-graining function is taken constant inside a 
volume V and zero outside it, provided that the interpar- 
ticlc separation is much smaller than the coarse-graining 
length scale (typically v^). 

The above expressions can be used to calculate the 
macroscopic fields from the microscopic ones (obtained 
e.g., from simulations or experiments), and compare them 
to the predictions of macroscopic models or direct exper- 
imental results. In order to close the set of continuum 
equations [Eqs. (3)] the stress and energy flux (the latter 
is not related to the considerations below) need to be ex- 
pressed as functionals of the pertinent macroscopic fields. 
As mentioned, such constitutive relations are often ob- 
tained empirically or conjectured. In some cases they are 
derived from the microscopic dynamics. The above exact 
expression for the stress field provides a framework for 
a systematic derivation of constitutive relations (as sug- 
gested for elastic networks in [27]). 

Here, wc arc concerned with the interpretation of ex- 
perimental data in terms of microscopic variables and macro- 
scopic fields. The fact that the contact stress includes a 
Slim over all contacts for each particle (i.e., even for very 
small CG scales the stress components correspond to spe- 
cific "averages" over the forces on each particle) already 
suggests that a "picture" of the forces in the packing does 
not correspond directly to the macroscopic stress field 
(they are certainly related, i.e., one would usually expect 
large stress components in regions where large force mag- 
nitudes are observed). In particular, as shown in Sec. 2, 
and further discussed below, force chains do not necessar- 
ily indicate macroscopic anisotropy or inhomogeneity. 
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Numerical results for model frictionless systems 

Consider a two-dimensional system of uniform disks (ar- 
ranged on a triangular lattice) subject to a vertical exter- 
nal force at the center of the top layer [26] . Experiments 
on such systems are described in [22,23]. Consider first 
the case of nearest neighbor harmonic interactions, i.e., 
the disks are coupled by equal linear springs (whose rest 
length is the diameter of a disk). Clearly, real cohesion- 
less particles do not experience any significant attractive 
interactions; however, there are a few insights to be ob- 
tained from the study of this system. Fig. 5 presents the 
forces in the system. Force chains arc evident. 

A contour plot of the "vertical stress component" azz 
[computed using Eq. (4)] for the same system is shown in 
Fig. 6 (with ^(r) = ■:;^e~^^'^^^'^^ , and w = d, the parti- 
cle diameter, i.e., a fine resolution). The force chains are 
not evident any more. The model described above cor- 
responds, in the continuum (long-wavelength) limit, to 
an isotropic 2D elastic medium [56]. The observed force 
chains, which break isotropy, can be attributed to the fact 
that the local environment of a particle in contact with a 
finite number of other particles cannot be isotropic. Under 
homogeneous macroscopic deformation, all forces would 
be equal in a lattice configuration. However, the concen- 
trated applied force yields an inhomogeneous deforma- 
tion, which leads to the local anisotropy being reflected 
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Fig. 5. Force chains in a 2D triangular lattice. A vertical force 
is applied at the center of the top layer. Line widths are pro- 
portional to the force magnitudes. Only the central part of the 
system is shown; reproduced from [26]. 




Fig. 6. Contour plot of hazz, corresponding to Fig. 5 {h is the 
slab height); reproduced from [26]. 



in the distribution of the forces. The elastic continuum 
description of the stress (to linear order in the strain) is 
isotropic, and cannot be expected to reflect this micro- 
scopic anisotropy. For small system sizes (in which the 
strain gradients on a particle scale are relatively large), 
this anisotropy can be observed in the stress field (a very 
clear example is shown in [26] for a macroscopically isotropic 
3D system, whose microscopic symmetry is cubic). These 
results, as well as those presented in Sec. 1 for disordered 
elastic systems, show that force chains do not necessarily 
indicate anisotropy or inhomogeneity of the material on 
sufficiently large scales; more importantly their existence 
does not require a non-elastic (microscopic) interaction. 

Note that only the forces between the particles and 
the floor (a single such force per particle) arc used in the 
calculation of the stress at the bottom of the packing. 
Hence on the bottom (but not in the bulk of the system), 
the spatial distribution of is equivalent (up to coarse- 
graining) to that of the microscopic forces. For sufficiently 
large systems, the distribution of forces on the bottom 
corresponds closely to the stress calculated using linear 
elasticity [26], even "almost without coarse-graining", i.e., 
for a microscopic CG scale. 



Fig. 7. Force chains in a 2D triangular lattice of 'one-sided' 
springs. A gravitational force has been applied in order to sta- 
bilize the system (the applied force is 150 times the particle 
weight); reproduced from [26]. 



A more realistic force model consists of 'one-sided' 
springs, i.e., springs that snap when in tension. Fig. 7 
presents the forces obtained for the same system presented 
in Fig. 7, but with 'one-sided' springs. Compared to the 
system of regular springs, the apphcation of the concen- 
trated force at the top of the packing leads to rearrange- 
ments in the contact network: some horizontal springs in 
the region under the point where the force is applied are 
disconnected (as also observed in [57] for a pile geome- 
try) but the force chains in both systems are qualitatively 
similar. The force distribution vs. the horizontal coordi- 
nate at different depths is in good agreement [26] with 
experiment [22,23]. For slightly disordered systems [26], 
the force chains are qualitatively similar, though some- 
what shorter. 

The corresponding vertical stress field u^z is shown in 
Fig. 8. The stress field in this case is clearly quite different 
from that obtained using 'two-sided', harmonic springs: 
the response for 'one-sided' springs is double-peaked. This 
is obviously related to the disconnected springs below the 
point of application of the external force. In [26], it has 
been shown that a model with harmonic springs in which 
the spring constant for the horizontal springs, is differ- 
ent from that of the oblique springs, corresponds (in 
the continuum limit) to an anisotropic clastic system. For 
sufficiently large K-ijKx, the response of such an elastic 
system has two peaks [26] (see also [18] for a more de- 
tailed analysis of the case of an infinite half-plane; the re- 
sults presented in [26] are for a finite slab on a rigid floor). 
The absence of horizontal springs corresponds to the limit 
K2/K1 00, the extreme anisotropic limit, which corre- 
sponds to an isostatic system. Note that the stress field, 
but not the displacement, depends only on K2/K1. In 
the case considered here, Ki = and K2 is finite; for 

— > 00, the rigid limit, the displacement is zero. The 
double peaked stress distributions are similar to those ob- 
tained from hyperbolic models. It follows that hyperbolic- 
like behavior can be obtained using an anisotropic (yet, 
still elliptic) elastic model (which becomes formally 'hy- 
perbolic' in the limit of very large anisotropy; see also [18, 
58]). 



Fig. 8. Same as Fig. 6, for the case of 'one-sided' springs (for 
which the forces axe shown in Fig. 7); reproduced from [26]. 



The 'stress-induced anisotropy' [26] observed in the 
case of 'one-sided' springs can be thought of as a non- 
linear extension of the linear clastic continuum behavior 
obtained in a network of harmonic springs. While the (pos- 
sibly position-dependent) elastic moduli in linear elastic- 
ity are time-independent material properties, a possible 
extension would be to introduce a stress history depen- 
dence of the elastic moduli (i.e., the anisotropy induced by 
the breaking of contacts in certain regions may be consid- 
ered a result of a tensile stress in those regions). A similar 
type of stress-induced anisotropy has been suggested in 
the context of plastic models for soil mechanics [59] . If the 
particle positions do not change significantly, so that only 
the contact network is modified in response to the applied 
stress, the behavior can possibly modeled as 'incremen- 
tally elastic'. Under certain condition (corresponding to 
plastic yield), the system would no longer be able to sup- 
port the applied stress without a major rearrangement of 
the particles. Incipient plastic yield may possibly be re- 
lated to a local extreme anisotropy typical of a marginally 
stable isostatic configuration. 

The force chains obtained both for the harmonic case 
and the 'one-sided' case are quite similar [26] to those 
observed experimentally [22,23], as an average over dif- 
ferent realizations. This averaging is required due to ex- 
perimental variations: although the particles are arranged 
on lattice, there is still some disorder present due to some 
variability in particle diameter, and possibly also in the 
contact properties [60]. Indeed, while perfect atomic lat- 
tices may be obtained at low enough temperatures, since 
all atoms of the same isotope are exactly identical, macro- 
scopic particles are never trolly identical, so that perfect 
periodicity can never be obtained. It also appears that the 
force chains obtained using the two models are quite simi- 
lar (similar chains arc also obtained in slightly disordered 
systems [26]). The stress field appears to be more sensitive 
to the anisotropy induced by the applied force. 

In the experiments reported in [20,21,24], the forces 
on the floor were measured. In [20,21], the width of the 
pressure probe (wliic;li would correspond to the CG scale 
of the measured stress) was 10 — 30 particle diameters. 
The bottom stress profiles measured are quite consistent 
with continuum elasticity (note that the depths of the 
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systems studied were 20 — 300 particle diameters). Ex- 
perimental deviations from the predictions of isotropic 
elasticity [21] can be reproduced by anisotropic elastic- 
ity [26]. Narrower or wider peaks than those obtained for 
isotropic systems can be obtained for small anisotropy, 
while for very large anisotropy, two peaks are expected 
(see also [18]). An additional possible cause for deviations 
from the isotropic elastic calculations presented in [21] is 
finite rigidity of the floor [56]. The more shallow systems 
used in the experiments may even be small enough for 
the finite size effects [26] to be significant. Any anisotropy 
in these experiments is obviously much weaker than the 
strong anisotropy observed in the model ordered system of 
'one-sided' springs. Several effects may explain this: first, 
the systems used in the experiment are highly disordered, 
so that inhomogeneous, random anisotropy may be ex- 
pected on intermediate (already macroscopic) scales, pre- 
sumably averaging out to an isotropic, or nearly isotropic, 
behavior at sufficiently large scales. In this case, the large- 
scale effect of contacts breaking due to applied forces would 
be significantly less pronounced than in the ordered sys- 
tem described above. A second possibility is the effect of 
frictional forces, which may either prevent contacts from 
breaking, or reduce the anisotropy of the response. Third, 
the model systems discussed above were unstressed before 
the application of the force, while the experimental ones 
are pre-stressed by gravity, which may compress some of 
the contacts such that the tension due to the applied force 
is insufficient to break them. 

In [24], individual forces on the floor were measured, 
and the results were averaged over realizations (which, as 
mentioned, is not necessarily equivalent to spatial coarse- 
graining). The regular packings used in [24] (FCC and 
HCP) are macroscopically anisotropic. The fact that some 
of the horizontal contacts (contacts among particles in the 
same layer) may be absent increases the anisotropy further 
(possibly in an inhomogeneous way; as mentioned above, 
a granular packing cannot be perfectly periodic). The ex- 
treme limit in which there are no such horizontal contacts 
corresponds to an isostatic system. Such anisotropy (pos- 
sibly further enhanced by the applied force) may explain 
the discrete peaks observed for relatively shallow systems 
composed of 9 layers of particles (and the fact that they 
appear to be consistent with a picture of "force propa- 
gation" appropriate for isostatic systems). However, for 
deeper system (about 20 particle diameters), there ap- 
pears to be a crossover to a smoother behavior, which 
should correspond to the crossover to the continuum limit 
(note that the depth of the systems used in [24] was smaller 
than the depth required in our calculations on 3D sys- 
tems [26] for reaching the continuum limit, so deeper sys- 
tems may still show dependence on the depth). 
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Effects of Friction 

As shown, some features of granular response may be re- 
produced using models employing frictionless and even 
harmonically interacting particles. However, it is clear that 
friction is consequential for granular materials. 



For frictional spheres, the microscopic description, as 
described in Sec. 2, must be extended to include (at least) 
the orientations of the particles, and interparticle torques 
in addition to the forces. The description of static and ki- 
netic friction requires the use of more complicated force 
models, which depend on the particle orientations and 
their relative tangential velocities, and possibly on the his- 
tory of contact deformation (see e.g., [37-39,61]). 

In experiments performed on regular 2D packings of 
photoelastic disks [23], the directions and "strengths" of 
the force chains observed upon application of a localized 
force to the top of the packing appear to depend quite 
strongly on the angle of the applied force with respect to 
the horizontal (in the following, all angles are given with 
respect to the horizontal). A particularly intriguing effect 
is that for some angles, force chains appear not only in 
the lattice directions (0, ±60°, ±120°, 180° for a triangular 
lattice), but also, apparently, in new directions which can 
be identified as ±30° (in fact, in individual realizations, 
rather than their average as reported in [23], it appears 
that force chains appear also at ±90°, i.e., the vertical di- 
rection [60]). These directions correspond to next- nearest 
neighbor directions in the triangular lattice. Since inter- 
actions among the particles only exist for particles in con- 
tact, there is no direct next-nearest neighbor interaction. 
The fact that the forces themselves, and not just the con- 
tact points, appear to be aligned with these ±30° direc- 
tions, suggests that frictional forces among the particles 
(tangential to the contact normals, which result in inter- 
particle torques) are necessary for obtaining forces (and 
chains) at angles different from the lattice directions. For 
an applied force at ±90°, it appears that the frictional 
forces are small enough such that the results obtained in 
this case [22, 23] are described quite well by a model of fric- 
tionless particles with linear force-displacement laws [26]. 

In order to elucidate the role of frictional forces and 

torques in the quasi-static response of granular materials 
in general, and in particular in order to gain an under- 
standing of the experimental results mentioned above [23] , 
we performed discrete element simiilations with normal 
and tangential linear spring-dashpot forces among the par- 
ticles (see e.g., [10,32]), possibly the simplest model for 
frictional disks. The simulation parameters were chosen to 
correspond to those of the experimental system [60]. Ex- 
perimentally, the force-displacement law for the photoe- 
lastic disks was found to be fit quite well by / oc ^^^^ [60] , 
as predicted by the standard Hertz theory for elastic ellip- 
soids in contact (see e.g. [36]), rather than the linear rela- 
tion (with logarithmic correction) expected for cylinders 
in contact [62], which appears to imply that the contact 
region between the "disks" is elliptic rather than rectan- 
gular. The simulation model described above employs a 
linear force-displacement law, so that an effective mean 
spring constant was estimated on the basis of the range 
of forces used in the experiments. The tangential spring 
constant was taken to be one-half the normal spring con- 
stant (a rough estimate consistent with the Hertz-Mindlin 
model [63] for oblique contact forces). The normal and 
tangential spring constants used are A:„ = BOOOfhg/R and 
kt = 1500fhg/R, where R and m are the mean particle 
radius and mass, respectively, and g is the gravitational 
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45° 30° 



Fig. 9. Force chains in 2D packings of slightly polydisperse frictional particles. A force F with magnitude 150 times the mean 
particle weight is applied to the particle at the center of the top layer. The angle of the force with respect to the horizontal is 
indicated below each picture. The same realization of the packing was used in all cases. The region shown is the central third 
of the upper half of the system. 



acceleration. The friction coefficient used is /it = 0.94 for 
particle-particle contacts and /z^*'' = 0.35 for particle- 
wall contacts. The systems studied here are composed of 
polydisperse disks, with radii distributed uniformly in the 
interval [R - SR, R], where 5R/R = 8 • 10~^ (i.e., a small 
polydispersity) . 

The system is first relaxed to a static state under grav- 
ity (until the total kinetic energy per particle is less than 
10~^mgR), and then relaxed again with an external force 
applied at the center of the top layer (in some cases the 
force was increased linearly with time from zero to prevent 
the "buckling" of the top layer which leads to major rear- 
rangements; these are beyond the nearly elastic behavior 
considered here). 

For comparison with the experiments presented in [23] , 
an external force of magnitude F = ISOmg was applied to 
the center top particle at angles of 15°, 30°, 45°, 60°, 75°, 
and 90°. The simulated systems consisted of 29 rows of 
80 particles, which is similar to the size of the systems 
used in the experiments [22,23,60]. Fig. 9 presents the 
forces obtained for different applied force angles. The same 
particle configuration was used in all cases. No significant 
particle rotation occurred except for the particles adjacent 
to the one on which the force is applied. For a force at an 
angle of 15°, buckling occurred in the top row, causing 
major rearrangements. Such buckling was also observed 
in the experiments, where it was apparently stabilized, 
limiting the rearrangements to a small region near the 
point of application of the force, but we have not been 
able to prevent major rearrangements in the simulation. 
As mentioned, tangential forces such as friction give rise to 
interparticle torques. Simulations with an applied torque 
(in addition to the applied force) show that this torque 
does influence the observed force chains [64]. 

The results are quite similar, qualitatively, to those ob- 
served in the experiment [23] . Note that the results shown 
in Fig. 9 are for a single configuration, while the results 
presented in [23] are for an average over configurations. 
The results obtained in simulations for different realiza- 



tions of the disorder are qualitatively similar [64]. The 
agreement of the results obtained using a relatively sim- 
ple force model with the experiments is encouraging. A 
more detailed study of the effects of friction on the forces 
and the stress field will be presented elsewhere [64]. 

6 

Concluding remarks 

We have shown that the seemingly inconsistent results of 
different kinds of experiments studying the static response 
of granular packings to a concentrated force can all be 
understood within the same framework of an essentially 
elastic (elliptic) picture once the distinction between forces 
and stress is made and the possible consequences of small 
system size, as well as anisotropy, are taken into account. 
The effect of applied stresses on the contact network may 
be modeled as a nonlinear, incrementally elastic model 
(which may be further extended to describe yielding). 

Somewhat surprisingly, many aspects of the response 
of such systems can be understood using models of fric- 
tionless particles. However, some effects do require the in- 
troduction of friction, as in the example of the force chains 
obtained for oblique applied forces described in this paper. 
We note, however, that the model for the friction used in 
the simulations described in Sec. 5 consists of tangential 
springs (with the additional Coulomb condition). This in- 
dicates that even for static frictional systems (below yield) 
an elastic continuum model, which probably includes ro- 
tational degrees of freedom (e.g., a Cosserat continuum 
model [65]), may be appropriate. 

Another important issue which requires further study 
is the effect of disorder, and the relation of spatial aver- 
aging to averaging over the disorder. 
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